#include <fstream>
#include <iostream>
#include <iomanip>
#include "TFile.h"
#include "TTree.h"
#include "TLine.h"
#include "TSystem.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TH1.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TArc.h"
#include "TEllipse.h"
#include "TClonesArray.h"

using namespace std;

void Ar(){

  int NO = 10;

  int ret = gSystem->Load( "../Objects/libObjectClass.so" );
  if( ret!=0 ){
    cout << "Failed to Load Library: ret = " << ret << endl;
    return;
  }

  TFile *fread=new TFile("/scratch/shinichi/OutputFiles/Ar_100keV_100ev.root");
  if(!fread){ cout<<"No Such File"<<endl; return;}

  double ThreHit = 5.;

  TTree *tree=(TTree*)fread->Get("tree");
  
  TClonesArray *JPMTa=new TClonesArray("PMT");
  TClonesArray *MPMT1a=new TClonesArray("PMT");
  TClonesArray *MPMT2a=new TClonesArray("PMT");
  TClonesArray *MPMT3a=new TClonesArray("PMT");
  TClonesArray *MPMT4a=new TClonesArray("PMT");
  TClonesArray *Particlea=new TClonesArray("Particle");
  TClonesArray *Neutrona=new TClonesArray("Particle");
  TClonesArray *Argona=new TClonesArray("Particle");
  TClonesArray *Gammaa=new TClonesArray("Particle");
  int nPhoton,nElectron;
  int nJhit,nM1hit,nM2hit,nM3hit,nM4hit;

  tree->GetBranch("JPMTb")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT1b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT2b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT3b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT4b")->SetAutoDelete(kFALSE);
  tree->GetBranch("Particleb")->SetAutoDelete(kFALSE);
  tree->SetBranchAddress("JPMTb",&JPMTa);
  tree->SetBranchAddress("MPMT1b",&MPMT1a);
  tree->SetBranchAddress("MPMT2b",&MPMT2a);
  tree->SetBranchAddress("MPMT3b",&MPMT3a);
  tree->SetBranchAddress("MPMT4b",&MPMT4a);
  tree->SetBranchAddress("Particleb",&Particlea);
  tree->SetBranchAddress("Neutronb",&Neutrona);
  tree->SetBranchAddress("Argonb",&Argona);
  tree->SetBranchAddress("Gammab",&Gammaa);
  tree->SetBranchAddress("nPhoton",&nPhoton);
  tree->SetBranchAddress("nElectron",&nElectron);
  tree->SetBranchAddress("nJPMThits",&nJhit);
  tree->SetBranchAddress("nMPMT1hits",&nM1hit);
  tree->SetBranchAddress("nMPMT2hits",&nM2hit);
  tree->SetBranchAddress("nMPMT3hits",&nM3hit);
  tree->SetBranchAddress("nMPMT4hits",&nM4hit);

  int nev = tree->GetEntries();
  cout<<"NEV = "<<nev<<endl;

  TGraph *ArgonXY[32];
  TGraph *ArgonXZ[32];
  TGraph *GammaXY[32];
  TGraph *GammaXZ[32];
  for(int i=0;i<32;i++){
    ArgonXY[i]=new TGraph();;
    ArgonXZ[i]=new TGraph();;
    GammaXY[i]=new TGraph();
    GammaXZ[i]=new TGraph();
  }

  TH1F *hNPhoton[8];
  TH1F *hNElectron[8];
  for(int i=0;i<8;i++){
    hNPhoton[i]=new TH1F(Form("hPhoton%d",i),Form("hPhoton%d",i),1000,1800,2800);
    hNPhoton[i]->SetXTitle("# of Photon");
    hNPhoton[i]->SetYTitle("NEV");
    hNElectron[i]=new TH1F(Form("hElectron%d",i),Form("hElectron%d",i),500,0,50000);
    hNElectron[i]->SetXTitle("# of Photon");
    hNElectron[i]->SetYTitle("NEV");
  }
  TH2F *hJPMTxy=new TH2F("hJPMTxy","hJPMTxy",150,-150,150,150,-150,150);
  hJPMTxy-> SetTitle("PMT Hit Distribution");
  hJPMTxy->SetXTitle("X (mm)");
  hJPMTxy->SetYTitle("Y (mm)");
  TH2F *hMPMT1xy=new TH2F("hMPMT1xy","hMPMT1xy",150,-150,150,150,-150,150);
  TH2F *hMPMT2xy=new TH2F("hMPMT2xy","hMPMT2xy",150,-150,150,150,-150,150);
  TH2F *hMPMT3xy=new TH2F("hMPMT3xy","hMPMT3xy",150,-150,150,150,-150,150);
  TH2F *hMPMT4xy=new TH2F("hMPMT4xy","hMPMT4xy",150,-150,150,150,-150,150);

  TH2F   *hPMTeff=new TH2F("hPMTeff","hPMTeff",30,-150,150,35,-650,-300);
  hPMTeff-> SetTitle("Combined");
  hPMTeff->SetXTitle("R (mm)");
  hPMTeff->SetYTitle("Z (mm)");
  TH2F  *hJPMTeff=new TH2F("hJPMTeff","hJPMTeff",30,-150,150,35,-650,-300);
  hJPMTeff-> SetTitle("J-PMT");
  hJPMTeff->SetXTitle("R (mm)");
  hJPMTeff->SetYTitle("Z (mm)");
  TH2F *hM1PMTeff=new TH2F("hM1PMTeff","hM1PMTeff",30,-150,150,35,-650,-300);
  hM1PMTeff-> SetTitle("Meg-PMT1");
  hM1PMTeff->SetXTitle("R (mm)");
  hM1PMTeff->SetYTitle("Z (mm)");
  TH2F *hM2PMTeff=new TH2F("hM2PMTeff","hM2PMTeff",30,-150,150,35,-650,-300);
  hM2PMTeff-> SetTitle("Meg-PMT2");
  hM2PMTeff->SetXTitle("R (mm)");
  hM2PMTeff->SetYTitle("Z (mm)");
  TH2F *hM3PMTeff=new TH2F("hM3PMTeff","hM3PMTeff",30,-150,150,35,-650,-300);
  hM3PMTeff-> SetTitle("Meg-PMT3");
  hM3PMTeff->SetXTitle("R (mm)");
  hM3PMTeff->SetYTitle("Z (mm)");
  TH2F *hM4PMTeff=new TH2F("hM4PMTeff","hM4PMTeff",30,-150,150,35,-650,-300);
  hM4PMTeff-> SetTitle("Meg-PMT4");
  hM4PMTeff->SetXTitle("R (mm)");
  hM4PMTeff->SetYTitle("Z (mm)");

  TH1F  *hJPMTdet=new TH1F( "hJPMTdet", "hJPMTdet",10,0,10);
  TH1F *hMPMT1det=new TH1F("hMPMT1det","hMPMT1det",10,0,10);
  TH1F *hMPMT2det=new TH1F("hMPMT2det","hMPMT2det",10,0,10);
  TH1F *hMPMT3det=new TH1F("hMPMT3det","hMPMT3det",10,0,10);
  TH1F *hMPMT4det=new TH1F("hMPMT4det","hMPMT4det",10,0,10);

  int NARGON,NGAMMA;
  int Jsum = 0; 
  int M1sum = 0; 
  int M2sum = 0; 
  int M3sum = 0; 
  int M4sum = 0; 
  int counterJ  = 0;
  int counterM1 = 0;
  int counterM2 = 0;
  int counterM3 = 0;
  int counterM4 = 0;

  double tQE = 0.025;

  for(int iev=0;iev<nev;iev++){
    JPMTa->Clear();
    MPMT1a->Clear();
    MPMT2a->Clear();
    MPMT3a->Clear();
    MPMT4a->Clear();
    Particlea->Clear();
    Neutrona->Clear();
    Argona->Clear();
    Gammaa->Clear();

    tree->GetEntry(iev);

    //----- NPhoton
    hNPhoton[0]->Fill(nPhoton);

    //----- PMT efficiency
    int ii = iev / 100;
    int iz = ii / 21 % 18;
    int ir = ii % 21;
    float zpos = iz*10.-545;
    float rpos = ir*10.-120;
    hPMTeff->Fill( rpos,zpos,(nJhit+nM1hit+nM2hit+nM3hit+nM4hit)/100.*tQE);
     hJPMTeff->Fill( rpos,zpos, nJhit/100.*tQE );
    hM1PMTeff->Fill( rpos,zpos,nM1hit/100.*tQE );
    hM2PMTeff->Fill( rpos,zpos,nM2hit/100.*tQE );
    hM3PMTeff->Fill( rpos,zpos,nM3hit/100.*tQE );
    hM4PMTeff->Fill( rpos,zpos,nM4hit/100.*tQE );

    Jsum  += nJhit;
    M1sum += nM1hit;
    M2sum += nM2hit;
    M3sum += nM3hit;
    M4sum += nM4hit;
    if( iev%100==99 ){
      hJPMTdet ->Fill(  Jsum/100.*tQE );
      hMPMT1det->Fill( M1sum/100.*tQE );
      hMPMT2det->Fill( M2sum/100.*tQE );
      hMPMT3det->Fill( M3sum/100.*tQE );
      hMPMT4det->Fill( M4sum/100.*tQE );
      if(  Jsum/100.*tQE>=ThreHit )  counterJ++;
      if( M1sum/100.*tQE>=ThreHit ) counterM1++;
      if( M2sum/100.*tQE>=ThreHit ) counterM2++;
      if( M3sum/100.*tQE>=ThreHit ) counterM3++;
      if( M4sum/100.*tQE>=ThreHit ) counterM4++;

      Jsum = 0; 
      M1sum = 0; 
      M2sum = 0; 
      M3sum = 0; 
      M4sum = 0; 
    }


    //----- PMT Hit
    for(int iJhit=0;iJhit<nJhit;iJhit++){
      PMT *jpmt = (PMT*)JPMTa->At(iJhit);
      hJPMTxy->Fill( jpmt->xpos,jpmt->ypos );
    }
    for(int iM1hit=0;iM1hit<nM1hit;iM1hit++){
      PMT *mpmt1 = (PMT*)MPMT1a->At(iM1hit);
      hMPMT1xy->Fill( mpmt1->xpos,mpmt1->ypos );
    }
    for(int iM2hit=0;iM2hit<nM2hit;iM2hit++){
      PMT *mpmt2 = (PMT*)MPMT2a->At(iM2hit);
      hMPMT2xy->Fill( mpmt2->xpos,mpmt2->ypos );
    }
    for(int iM3hit=0;iM3hit<nM3hit;iM3hit++){
      PMT *mpmt3 = (PMT*)MPMT3a->At(iM3hit);
      hMPMT3xy->Fill( mpmt3->xpos,mpmt3->ypos );
    }
    for(int iM4hit=0;iM4hit<nM4hit;iM4hit++){
      PMT *mpmt4 = (PMT*)MPMT4a->At(iM4hit);
      hMPMT4xy->Fill( mpmt4->xpos,mpmt4->ypos );
    }
    


  }

  //----- Draw Part

  TLine *lFid1=new TLine(-133,-360,-133,-550);
  TLine *lFid2=new TLine(-133,-550,  94,-550);
  TLine *lFid3=new TLine(  94,-550,  94,-360);
  TLine *lFid4=new TLine(  94,-360,-133,-360);

  float max1 = 5.;
  TCanvas *c1=new TCanvas("c1","c1",800,600);
  hPMTeff->SetMinimum(0.);
  hPMTeff->SetMaximum(max1);
  hPMTeff->Draw("colz");
  lFid1->Draw();
  lFid2->Draw();
  lFid3->Draw();
  lFid4->Draw();

  //___ Photon Distribution
  float max2=2000;  
  TCanvas *c2=new TCanvas("c2","c2",680,600);
  hJPMTxy->SetMinimum(0.);
  hJPMTxy->SetMaximum(max2);
  hJPMTxy->Draw("colz");
  hMPMT1xy->SetMinimum(0.);
  hMPMT1xy->SetMaximum(max2);
  hMPMT1xy->Draw("samecolz");
  hMPMT2xy->SetMinimum(0.);
  hMPMT2xy->SetMaximum(max2);
  hMPMT2xy->Draw("samecolz");
  hMPMT3xy->SetMinimum(0.);
  hMPMT3xy->SetMaximum(max2);
  hMPMT3xy->Draw("samecolz");
  hMPMT4xy->SetMinimum(0.);
  hMPMT4xy->SetMaximum(max2);
  hMPMT4xy->Draw("samecolz");
  TArc *aInner=new TArc(0.,0.,133.,45.,315.);
  TLine *lInner=new TLine(133./sqrt(2),-133./sqrt(2),133./sqrt(2),133./sqrt(2));
  aInner->SetFillStyle(0);
  aInner->SetNoEdges(1);
  aInner->Draw();
  lInner->Draw();

  //___ Photon eff
  gStyle->SetOptStat(000000);
  float Zmax = 10;

  TCanvas *c3=new TCanvas("c3","c3",1200,900);
  c3->Divide(3,3);
  c3->cd(5);
  hJPMTeff->SetMinimum(0.);
  hJPMTeff->SetMaximum(Zmax);
  hJPMTeff->Draw("colz");
  lFid1->Draw();
  lFid2->Draw();
  lFid3->Draw();
  lFid4->Draw();
  c3->cd(3);
  hM1PMTeff->SetMinimum(0.);
  hM1PMTeff->SetMaximum(Zmax);
  hM1PMTeff->Draw("colz");
  lFid1->Draw();
  lFid2->Draw();
  lFid3->Draw();
  lFid4->Draw();
  c3->cd(1);
  hM2PMTeff->SetMinimum(0.);
  hM2PMTeff->SetMaximum(Zmax);
  hM2PMTeff->Draw("colz");
  lFid1->Draw();
  lFid2->Draw();
  lFid3->Draw();
  lFid4->Draw();
  c3->cd(7);
  hM3PMTeff->SetMinimum(0.);
  hM3PMTeff->SetMaximum(Zmax);
  hM3PMTeff->Draw("colz");
  lFid1->Draw();
  lFid2->Draw();
  lFid3->Draw();
  lFid4->Draw();
  c3->cd(9);
  hM4PMTeff->SetMinimum(0.);
  hM4PMTeff->SetMaximum(Zmax);
  hM4PMTeff->Draw("colz");
  lFid1->Draw();
  lFid2->Draw();
  lFid3->Draw();
  lFid4->Draw();


  //___ hNphoton
  TCanvas *c4=new TCanvas("c4","c4",800,600);
  hNPhoton[0]->Draw();


  //___ 
  TCanvas *c5=new TCanvas("c5","c5",800,600);
  TH1 *frame5=gPad->DrawFrame(0,0,10,380);
  frame5->SetXTitle("# of Detected Photon");
  frame5->SetYTitle("NEV");
  hJPMTdet ->Draw("same");
  hMPMT1det->SetLineColor(2);
  hMPMT1det->Draw("same");
  hMPMT2det->SetLineColor(3);
  hMPMT2det->Draw("same");
  hMPMT3det->SetLineColor(4);
  hMPMT3det->Draw("same");
  hMPMT4det->SetLineColor(6);
  hMPMT4det->Draw("same");

  cout<<" Jhit>"<<ThreHit<<": "<<counterJ<<" / "<<378<<endl;
  cout<<"M1hit>"<<ThreHit<<": "<<counterM1<<" / "<<378<<endl;
  cout<<"M2hit>"<<ThreHit<<": "<<counterM2<<" / "<<378<<endl;
  cout<<"M3hit>"<<ThreHit<<": "<<counterM3<<" / "<<378<<endl;
  cout<<"M4hit>"<<ThreHit<<": "<<counterM4<<" / "<<378<<endl;

}
